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I. INTRODUCTION 



D. A. Quint 1 and J. M. Schwarz 1 
1 Physics Department, Syracuse University, Syracuse, NY 13244 
(Dated: February 23, 2011) 

Actin cytoskeletal protrusions in crawling cells, or lamellipodia, exhibit various morphological 
properties such as two characteristic peaks in the distribution of filament orientation with respect 
to the leading edge. To understand these properties, using the dendritic nucleation model as a basis 
for cytoskeletal restructuring, a kinetic-population model with orientational-dependent branching 
(birth) and capping (death) is constructed and analyzed. Optimizing for growth yields a relation 
between the branch angle and filament orientation that explains the two characteristic peaks. The 
model also exhibits a subdominant population that allows for more accurate modeling of recent 
measurements of filamentous actin density along the leading edge of lamellipodia in keratocytes. 
Finally, we explore the relationship between orientational and spatial organization of filamentous 
actin in lamellipodia and address recent observations of a prevalence of overlapping filaments to 
branched filaments — a finding that is claimed to be in contradiction with the dendritic nucleation 
. model. 

<D 
<N 

The process of cell motility involves a number of components — the actin cytoskclcton, the cellular membrane, an 
assortment of actin-binding proteins, molecular motors and intcgrins — that assist the cell in changing shape so that 
, it can move in a particular direction Naturally, one assumes that the interplay between the various components 
has been tuned to form structures that optimize for efficient motility. To test this assumption quantitatively is not 
necessarily an easy task given the dynamically complex structures that emerge as a cell crawls. However, theoretical 
descriptions of complex biological systems rooted in simplicity may helpto identify key interactions so that the 
I ' sophistication of cell motility may be better quantitatively understood [H,Q- 

In order for the cell to crawl in a particular direction, the cell extends itself. This extension, otherwise known as 
the lamellipodium, is facilitated by the growth and restructuring of the actin cytoskclcton. Over the past ten years or 
, so, the dendritic nucleation model has emerged as the dominant conceptual picture of this reorganization [HQ. The 
' dendritic nucleation model asserts that cytoskeletal growth is initiated by membrane-bound proteins, such as PIP2, 
t-H \ that activate WASP. WASP, in turn, activates Arp2/3, a protein that nucleates new filaments from preexisting ones. 
0^ • At the point of nucleation, the branch angle takes on a somewhat regular angle of 70 degrees with respect to the 
mother filament [1, H|. This nucleation takes place at/near the cell membrane and leads to a tree-like, or dendritic 
structuring of the actin cytoskclcton. 

New filament growth must be accompanied by some system of regulation, since unregulated birth of branched 
actin filaments can lead to a redundant use of finite resources in a cell. It has been shown in purified reconstituted 
systems that Arp2 /3 and G-actin alone are insufficient for motility 0, H[ . Additional regulation of the existing actin 
cytoskeleton is required for rapid G-actin turnover. This regulation is partially assisted by capping proteins, which 
attach to the plus ends of filaments and stop polymerization. In other words, the filament dies. Also, filaments further 
back from the leading edge debranch and get severed, eventually becoming part of the finite pool of G-actin. All of 
these processes are qualitatively described by the dendritic nucleation model. For a unified quantitative description 
d ' of such processes see Ref . [9| . 

Using the branching (birth) and capping (death) processes as a basis for formation of lamellipodia, we give quan- 
titative evidence to support the notion that form/morphology of the dendritic network is optimized to facilitate cell 
motility. More specifically, we propose a mean field model for the birth and death rates as a function of filament 
orientation with respect to the leading edge. We then optimize for filament reproduction at the leading edge, which 
provides an optimal relation between the branch angle and the angle with respect to the leading edge that agrees with 
experimental observation 0] . 

We must point out that there exists an earlier mean field model, which predicts the same optimal relation between 
the branch angle and the angle of orientation with respect to the leading edge. However, the earlier model has a 
different physical basis [loj . In keeping with the scientific method, we study further implications of the two models 
in order to make other retrodictions/predictions to distinguish them. For example, by studying the implications of 
our model on the spatial organization of filaments, we propose a new shape for the filament density along the leading 
edge. The shape may account for an observed "excess" filament density along the outer edges of the lamellipodia 
beyond what current modeling predicts [ll| . We also make comparisons with a more recent orientational model [l2j ■ 
Finally, our analysis of spatial information allows us to investigate a recent experimental study of lamellipodia 
made by Urban and collaborators using electron tomographic images of cytoskeletal networks [l3[ . They found that 
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overlapping actin filaments were much more prevalent than branched filaments. Based on this data, they proposed an 
alternate model for the reshaping of actin filaments near the leading edge — that polymerization and cross-linking are 
the main ingredients for cellular extension and not Arp2/3, which is relegated to a non-branching nucleating agent of 
new filaments just as dimerization nucleates new filaments. We implement a discrete, spatial simulation of our model 
to measure, for example, the ratio of overlaps to branch points to determine if the prevalence of overlaps rules out the 
dendritic nucleation model. We also use the full two-dimensional spatial information of the filament tips to discuss 
implications for the buckling of the network. 

The paper is organized as follows: Section II introduces and analyzes the mean field, orientational birth-death 
model. Comparisons with earlier mean field models, as well other generalizations, are addressed. Fluctuations about 
the mean field solutions are investigated. Section III studies the coupling between the orientational degrees of freedom 
and the spatial degrees of freedom, such as analyzing the full two-dimensional information of filament positions via 
discrete simulations. Section IV discusses the implications of our results. 



II. MEAN FIELD MODELS 



A. Collision-based model 



Actin filaments contain an inherent polarity where actin monomers associate with the plus end of the filament and 
dissociate from the minus end This polarity allows for directed assembly such that the cell can extend itself 

in a particular direction. While extension via polymerization is one mechanism for extension, the dendritic nucleation 
model [HH! asserts that extension via nucleation of branched filaments off pre-existing ones is also important. Support 
for the dendritic nucleation model has come about, for example, from electron micrograph images of branched 
actin networks in lamellipodia, from the knocking out of Arp2/3 preventing the formation of lamellipodia [17| and 
from reconstituted systems of purified proteins 0. In these reconstituted systems, motility can be induced by using 
a small number of purified proteins combined in vitro which can reach speeds, for optimal concentrations, of 2.2 \im 
minT 1 . It was observed that motility cannot occur with activated Arp2/3 alone. Additional proteins, which facilitate 
a steady state of G-actin concentration, are essential Q- These proteins are capping proteins, which cap polymerizing 
plus ends, and ADF/Cofilin which cuts actin filaments. Both proteins help replenish the G-actin pool. 

To test some of the dendritic nucleation model assertions, let us construct a mean field model with branching and 
capping and investigate various experimental consequences. As for the branching, in vitro studies suggest a preferred 
angle of 70° with respect to the plus end of the mother filament (if. Therefore, for now, we assume that the branching 
angle is some fixed angle ip from the plus end. We will also assume that Arp2/3 branches off the side of pre-existing 
filaments with a preference towards the plus end as has been observed experimentally [l8j . Moreover, we assume that 
the nucleation of a branch occurs at and/or very near the membrane. Of course, if branching takes place only at 
the membrane, then the initial structure of the network is dictated by the shape of the membrane. If the Arp2/3 is 
released from the membrane upon activation and then collides (binds) with actin filaments, then the network structure 
is less dependent on the shape of the membrane. Recent experiments observed space-filling polymerization of filopodia 
into gaps between the edge of the network and the membrane fl9j . Such an observation has to yet found with Arp2/3 
nucleation, however. 

So, assuming that side-branching occurs and that the branch is nucleated at/near the membrane, the branching 
probability depends on the orientation of the filament. The more the filament is parallel with the leading edge, 
the higher the cross-section for collision between the globular Arp2/3 and the one-dimensional filament and, hence, 
nucleation. More precisely, the branching rate contains a | sin((9)| dependence, where 9 = 0° is normal to the leading 
edge. See Figure 1. 

As for the death rate, filament plus ends get capped at a rate c. We will not assume any angular dependence 
for the capping rate. The capping protein-plus end binding is a globular-to-globular collision. Furthermore, the 
task of the capping protein is to regulate the length of filaments such that growth is channelled into developing new 
filaments and not into extending pre-existing ones [2(|. Elongating pre-existing filaments leads to a system longer 



filaments on average, which are more susceptible to buckling [2 lj . Channeling new filament growth should, therefore, 



be independent of filament orientation. In addition, channelling growth into branches allow for further spreading of 
lamellipodia, which increases cell contact with the surface in order to build more focal adhesions. 

Based on the above assumptions, we construct a set of kinetic equations that take into account the orientation of 
filaments, which branch off prexisting filaments and get capped. We restrict ourselves to —90° < 9 < 90° since we 
are only interested in "forward" growth. We first consider ip > 45° and 0° < 9 < ip. In this regime, there are are 
two populations of filaments, filaments oriented at angle 9, denoted by n\, and filaments oriented at an angle 9 — "0, 
denoted by n-i- (There is a reflection symmetry about 9 = 0°. We will only deal with 0° < 9 < 90° and use the 
reflection symmetry to extend our results to —90° < 9 < 0°.) The kinetic equations for this first case are 
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FIG. 1: Schematic of the orientation of a branched filament (B) in relation to its mother filament (M) and the horizontal line 
denotes the leading edge. The dashed line represents the normal to the leading edge. 



= !lsin(0- VOK -cm (1) 

and 

^ = ^|sin(0)|ni -cn 2 , (2) 

where b denotes the magnitude of the branching rate. The factor of 1/2 is because branching on the "backside" of 
the mother filament is a less- likely collision given the activation of Arp2/3 at the membrane and we do not consider 
it here. 

We now rephrase famous "the form follows function" optimization guideline as a population biology problem (22"142~4| . 
We assume the cytoskeletal system is maximizing for "population" growth so that the cell can extend itself efficiently. 
To determine this maximal growth, we compute the eigenvalues for the above set of equations and determine the 
relationship between 9 and ip such that the largest of the two eigenvalues is maximized (and positive). The eigenvalues 
for the above set of equations are 



Ai, 2 = -c± -V| sin(0)|| sin(0 - ^)\. (3) 

It is easy to see that the largest eigenvalue is maximized when 9* = ip/2. Of course, 9* = —ip/2 is another optimal 
solution via symmetry. 

Next, we investigate ip < 9 < 90° (and xp > 45°). In this second regime, there are three orientations of filaments 
with 713 denoting filaments oriented at 9 — 2ip. The set of kinetic equations for this second case are 

= ^|sin(0-V>)K -cni, (4) 
^ = b -\sm(9)\n 1 + ^\ S m(9-2^)\n 3 -cn 2l (5) 

-7T = ^|sin(6»- V)|n 2 -cn 3 . (6) 

Once again, assuming activation of Arp2/3 at the membrane, nucleation is unlikely to occur on the "backside" of 
a mother filament with respect to the membrane and the kinetic equations become 

dm 

= ~c ni , (7) 

^ = ||sin(e)|ni + ||sin(ff-2^)|n 3 -cn 3> (8) 

dnz b. . , .. , . 

— = -|sin(0- ip)\n 2 -cn 3 . (9) 
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FIG. 2: Depiction of the optimal orientation (red) and the two suboptimal ones (blue, green). 



The rii population eventually dies off such the above equations simplify further to 



dn 2 

~~dt 

dn 3 

dt 




2ip)\n 3 - cn 2 , 



tf))\n 2 - cn 3 . 



(10) 



(11) 



With the transformation of 8' = 8 — ip, we map back to the first set of kinetic equations such that the largest, 
positive eigenvalue occurs again at 8'* — ±-0/2. If ip = 45°, there exists another optimum at 8* — 67.5°. However, 
this initial orientation will die away and the ±22.5° will survive. So we have the same optimization as in the first 
case, i.e. it is redundant. This result is different from the initial model where the optima occur at 8 = 0°,±ip. 
Of course, for 8 = 0°, the critical buckling load is the smallest, i.e. filaments are more susceptible to buckling. While 
this property is not optimal for rheology, bundled filaments can increase the critical buckling load [25j . For our model, 
there can be no optimum at 8 = 0°. 

If we increase ip beyond 60°, the second optimization peak is outside of the range of interest (ip < 8 < 90°). However, 
consider ip = 70°. As 8 increases from 70° to 90°, the reproductive growth enhances monotonically. For 8 = 90°, 
the daughter filaments are subsequently oriented at 20° and —50°. While these initial 90° filaments are not precisely 
optimized for growth, their growth is maximum for the interval, ip < 9 < 90°. Therefore, we deem these 20° and —50° 
filaments as suboptimal orientations. Figure 2 depicts the optimal orientations and the suboptimal, or subdominant, 
orientations. We conjecture that the subdominant orientations of filaments may serve as reinforcements for cross- 
linking. Depending on the initial spatial arrangement and orientation of the filaments, the subdominant orientations 
may help to increase overlaps and spreading. It is interesting to note that only for ip > 60°, the second redundant 
optimum is removed. The observed branch angle is reasonably close to this value. 

Ex per iments on keratocytes have measured the distribution of orientation of filaments normal to the leading 
edge [10(. There are two maxima in the distribution occurring at ±35°. Assuming that the branch angle is indeed 
70°, our optimization analysis provides an explanation for this experimental finding. However, we should mention 
that Koestler and collaborators have conducted a more recent experiment on the orientation of filaments (2(|. They 
observed a broad distribution between the angles of —75° to 75°. One could argue that the subdominant orientation 
of filaments could account for further spreading of the distribution and, therefore, perhaps the more recent data. 



How does the above model compare with the one constructed and analyzed by Maly and Borisy [10(? The 
Maly/Borisy model is consistent with the Brownian ratchet model [13, [H[ for filament elongation near a mem- 
brane. In the Brownian ratchet model, leading edge filaments polymerize only if there exists enough space between 
the membrane and the tip of the filament. As the filaments fluctuate, transient gaps open up between the filament 
and the membrane, allowing actin monomers to attach to the plus end of the filament. Once the filament bends back 
to its original straight configuration, it is now longer and, therefore, pushes against the membrane moving it forward. 
This process, however, is limited to the size of the fluctuations that occur between the membrane and the tip of the 
leading edge filament. In support of this notion, experiments involving changing the membrane tension have shown 
that there exists an inverse relationship between the lamellipodial extcntion velocity and the apparent membrane 
tension [29| . However, more recent experiments suggest more complicated mechanisms may be at play (30|. 



B. Comparison with other models 
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Given the space limitation between the membrane and the fluctuating tip, there is an orientational degree of freedom 
that the filaments can exploit in this polymerization process. By varying the angle at which the tip makes with the 
membrane initially, the amount of space between the two can change if the membrane moves forward over a time t at 
a velocity v mem . In particular, 5k p p m pcos(9) = v mem , where S is the G-actin diameter, k p is the polymerization rate, 
p m is the G-actin concentration, and p is the probability that the filament tip is not obstructed by the membrane. 
Maly and Borisy assert that the capping of a filament is only possible if the growing filament tip is not obstructed 
by the membrane, hence capping occurs at a rate cp. Since p oc 1/ cos(#), the larger 9 is, the more likely the filament 
will be capped. 

As for the branching, in the Maly/Borisy model, the branching rate does not contain any angular dependence. In 
other words, the kinetic equations read 

dni b cp 

~df = ^{f) ni (12) 

dn 2 b cp 

-JT = o n i Ta 7T cri 2, (13) 

at 2 cos(9 — ip) 

where po = v mem /6p m k p . Maly and Borisy show that for cos(?/>) < p$ < cos(V>/2), the optimal relation between 9 and 
ip is, as above, 9* = ±■0/2. However, for p < cos(-0), the optimal orientations are zero and ±ip. 

A more recent orientational model assumes a ^-independent, zcroth-order branching rate, a 9- independent, first- 
order capping rate, and a ^-dependent outgrowth rate that kills single filaments outgrowing the bulk of the network [l2T | . 
The model exhibits two different, stable patterns, the same two exhibited by Maly and Borisy [Ioj |. 9 = ±0/2 or 
9 = 0, ±"0. The two patterns cannot coexist. Parameters such as the capping rate determine which pattern prevails. 
The authors argued that their model can explain the experimentally observed load-dependence of the network velocity 
at a given force [3lj . Our model cannot exhibit the latter pattern and our subdominant pattern for ip > 60° does 
coexist with the primary, or optimal, one. 

C. Generalized Birth/Death rates 

While each mean field model has a different physical basis, the selection criterion for maximal growth yields the 
same optimal relationship, 9* = ±ip/2. How generic is this result? To begin to answer this, we consider the most 
general version of our population equations such that both the birth-rate and the death-rate depend on the orientation 
of the branched filaments. Therefore, we begin with 



dri\ 



fli(0,VOna+.Di(0,^)ni ( 14 ) 
B 2 (6,i>) ni +D 2 (6,ip)n2, (15) 



for 0° < 9 < 90°. We define the matrix, Q = Q(0, ip) such that we can represent the set of linear coupled equations 
vectorially as n = Qn, where 

'Di(M) Bx(M) 
Defining Q = Q/Det[Q], the eigenvalues of Q are given by 



Q - (b 2 (0,V) D 2 {9^))- (16) 



Tr[Q] ± ^/Tr[QP - 4 
*+,- n ■ ( 17 ) 



With this result, three scenarios emerge: 



Condition Eigenvalues 

(1) \Tr[Q}\<2 A'+,_^C 

(2) |Tr[Q]| = 2 (A+ = A-)^R (18) 

(3) \Tr[Q} \ > 2 (A+ > A_) -> R 
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To determine the largest, real eigenvalue, we focus on condition 3. Dropping the +,- notation, the optimization 
condition is determined by 

dgX = Tr[dQ] (l + Tr[Q] ) = (19) 

such that 

Tr[d e Q] - 0. (20) 

In other words, the optimization condition occurs when the matrix <9gQ is rendered traceless. One, of course, also 
needs to evaluate the second derivative to check for a maximum. 

With the help of the Jacobi formula for the derivative of the determinant of a matrix and using the linearity of the 
trace operator, the optimization condition for Q must satisfy (assuming the trace of <9#Q is zero), 

IMQ-^Q] = 0. (21) 

If we analyze the case where the two death rates are 0- independent, then the optimal condition is 



dg(B 1 (9,ip)B 2 (9,ip)) = 0. (22) 

For example, if Bi(9,ip) is V'-independent and B2(9,ip) = B\(9 — ip), then the optimal condition is 9* = ip/2 as 
long as B\{9) is an even function (provided the second derivative is negative at that point). If B\ is a trigonometric 
function, then the periodicity should not be too small such that other maxima appear within the 0° < < 90° range. 
So, B\ = cos(0) yields the same optimal relation between 8 and ifj as would many other functions. It is possible to 
broaden this analysis. We leave this for future work. Our point now is that the optimal finding of 9* = ±-0/2 alone 
does not necessarily justify the model. One needs to explore further implications of the model in order to distinguish 
it from other potential models. We shall pursue this tact in the next section. 



D. Fluctuations 



Is the optimal relationship between 9 and tp robust in the presence of fluctuations? To answer this question, following 
Maly and Borisy [3, we assume that the angle between the two types of filaments exhibits Gaussian fluctuations 
with a mean of ip and a variance of a 2 . If we define n(9,t) as the density of filaments at the leading edge at time t 
and orientation 9, then the dynamic equation for n(0, t) (for — ip < 9 < ip) is given by 

dn(9,t) 6|sin(0)| f' 4 ' , c'+^-<» 2 v' . 
at v 87t<7 J-,^ 

If we assume n(9, i) = N(t)q(6), we arrive at 

|sin(0)| , (.e'+i>-<» 2 (9' -^-0)2 . . . , , 

71^7 \.} e — ^~ + e — ^)^ e ) de = yd )* ( 24 ) 

where y is now an unknown eigenvalue such that the above assumption is justified. We use the quadrature method 
to numerically solve for q{9) for different values of a. Figure 3 depicts the results. The maximum of q(9) correspond 
well with the largest, positive eigenvalue found previously. As a increases, the maxima remain robust, but are become 
less pronounced. These results indicate that the optimal relation of 0* = ±-0/2 is robust to fluctuations. 

It is certainly worth comparing this result with the fluctuation results of Maly and Borisy [l(| ■ The steady state 
orientation in the presence of noise here is very similar to the Maly/Borisy model [10j | , at least for intermediate values 
of po- This new computational result is, therefore, somewhat nontrivial given that one would assume the fluctuations 
to be sensitive to the details of the underlying kinetics. Further investigation along the lines of Section lie is needed 
to pursue understanding of this possible gencricity despite differing details of the kinetics. 
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FIG. 3: Plot of q(0) for 4> = 70°. For a = 4°, ^ = 0.269, for cr = 8°, = 0.254, and for a = 12° , ^ = 0.213. 

III. ORIENTATION INFLUENCING SPATIAL ORGANIZATION 
A. Filament density profile along the leading edge 

Optimization for growth in lamellipodia leads to a relationship between the branch angle ip and the orientation 
of filaments relative to the leading edge, or 9. To date, there exist three models, each rooted their in own physical 
basis, that yield 0* = ±ip/2. In order to further differentiate between these models, we investigate the distribution of 
filament tips along the leading edge. 

Previous work investigating the filament density along the leading edge has invoked the following set of assump- 
tions [ill HH HH • Filaments are either oriented with +35° or —35° with respect to the leading edge. Their respective 
densities along the leading edge x are denoted by p + (x,t) and p~(x,t). These filaments undergo lateral flow in their 
respective directions. Filaments with either orientation can spawn filaments with the opposite orientation (± — > =p) 
from their own. Also, both types of filaments can get capped. Therefore, the equations for both filament densities 
along the leading edge, whose position is denoted by x, are 



with B = J 2 l dx(p + (x) + p~(x))i where L is the length of the lamcllipodium and v is the lateral flow speed, which 

2 

is proportional to the speed of the crawling cell. 

Previous analysis of the above equation yields a total filament density in steady state that is peaked at the center of 
the cell, provided the filament density at the edges is sufficiently small. More specifically, for the boundary conditions, 
p + (-L/2,t) = and p~(L/2,t) = 0, p+(x,t -> oo) + p-(x,t — > oo) = § cos(^). Therefore, near the center the 
profile is an inverted parabola (ill. l33l] . If the boundary condition is adjusted to a higher concentration, eventually, the 
inverted parabola becomes a parabola with the total filament density higher at the sides than in the center. Assuming 
j- << c, we define a dimensionless time, r = ct, a dimensionless position s = x/L, and dimensionless densities, 
p^ 1 = p ± ^-- l then Eq. 25 becomes 

dp^ 8 , V I . 1 

-h=^-L^ + B pT - p ■ (26) 

If the boundary conditions at the ends of the lateral extent demand a large enough density, then the system will not 
be able to sustain the peak in density at the center of the leading edge. The larger the branching rate, the higher the 
allowed density at the ends can be with the system still sustaining an inverted parabola. 

The inverted parabola in filament density along the leading edge is observed in experiments near its center fllj ) . 
However, there is an excess of the filament density towards the sides of the leading edge (-L/2 and L/2) that appears 
to be flat. This excess has not been accounted for in the current model. In light of the collision-based model introduced 
here, we propose that for ip < < 90° with ip ~ 70°, the subdominant pattern, whose growth is prevented from being 
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fully optimized so as not to form a redundant pattern, may account for this excess. The axis of the subdominant 
pattern is at 9 = —15° as opposed to 6 = 0°. There exists another pair centered at 9 = 15° by symmetry. 

To test this proposal, we take the simplest approach by constructing the following six equations taking into account 
the two center populations (as before) and the two respective pairs of subdominant, or "off-center" , populations. For 
now, each respective pair of populations is not coupled to any other respective pair. Each pair occupies it own region 
along the lateral extent of the lamellipodia. We denote pt as the original set, pt as those directed toward the left 
side of the leading edge (from the birdseye perspective of cell) , and pt as those directed toward the right side of the 
leading edge to arrive at 



dpt 
dt 






-cpf 


dpf 
dt 




' ~Bi 


-c P f 


dpf 
dt 


= T Tx {vp r )A 


b 2 T 


- cpf 



(27) 
(28) 
(29) 

— -I- — — 4- _ 

where B c = §\ dx(pj (x) + p~ (x)) , Bi = J_l dx(pj (x) + p ; - (x)), and B r = J Q 2 dx(pj (x) + p^ (x)). Note the different 
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spatial regions for each respective pair of populations. Of course, the delineation is not so clear cut in practice. Also, 
b 2 < b to allow for a slight decrease in branching at the edges of the leading edge. Finally, we assume v does not vary 
between the different pairs. 

To solve for the steady state filament density distribution, we use the following boundary conditions. We set 
p+(-L/4,t) = 0, p~(L/4,t) = 0, p+(0,i) = 0, p~(L/2,t) = p Q , p^(0,t) = 0, pf{-L/2,t) = p . As for the asymmetric 
boundary conditions on pt and pf , it is reasonable to assume that near the lateral center of the leading edge, the 
density of /9+(0, t) = 0. However, towards the sides of the leading edge, p~ (L/2, t) may not necessarily vanish as there 
may be some skewing of the axis along which the subdominant populations are propagating due to the focal adhesions. 
The same goes for pf '. For the symmetric boundary conditions for pt, the same cosine steady state solution exists 
as before (only over a smaller interval). For the pt and pf populations, the steady state solutions are sinusoidal. (If 
pt(0,t) = pi « 1, then the steady state solution is a linear combination of sine and cosine). These solutions are 
plotted in Figure 4. We use j- = 0.01, b 2 = 0.6b, and pa = 0.3. Typical lateral speeds are of order 0.1 microns/sec 
for fast-moving keratocytes, typical lengths of leading edges are tens of microns and typical capping rates are tenths 
per second to per second 0, HH . We have also checked these solutions numerically. 

When we sum up the various densities to arrive at the total density, we can account the observed excess filament 
density at the lateral edges of the lamellipodium while still having the overall proper shape found experimentally 
near the center of the system. Indeed, there is a small dip in the center of the system, this dip may be difficult to 
observe experimentally and is presumably washed out once noise and other details are incorporated into the modeling. 
For instance, the revised model may be further updated to included coupling between the different populations via 
branching in terms of the overlapping regions. We leave this for future work. We should also note that the two 
previous orientational models, at least for 9* = 0°,±ip, would yield a filament density distribution that is sensitive to 
the initial distribution of filaments along the leading edge since the 9 — 0° does not flow laterally. Such a sensitivity 
should be investigated in order to rule out the possibility of the 9 = 0°, ±t/j pattern. 



B. Two-dimensional, discrete simulation 



To further study how filament orientation affects the spatial distribution of filaments, we construct a two-dimensional 
kinetic simulation with explicit filaments. A two-dimensional approach is reasonable given that lamellipodia are 
typically flat structures with a thickness of approximately 100 nanometers, extending several microns into the body 
of the cell and approximately 10 microns across the cell. The simulation algorithm is as follows: 

(0) Initialization: A filament is initialized at the origin of the system with an angle 9 and a length 100 nm. 

(1) Branching: A random number, r, is chosen from the sine distribution. Should r < sin(#) (with 9 < 90°), a branch 
point is chosen along the initial filament. Where the branch point occurs is uniformly chosen over some part of the 
current filament length as measured from the growing end. The length is denoted by /. This constraint restricts the 
branching to occur near the leading edge of the network. The branch filament emerges at an angle ip with respect to 
the mother filament. Gaussian fluctuations about the branch angle, with variance a 2 , are also studied. 

(3) Capping: A random number, s, is chosen uniformly between zero and unity. If s < c, the filament gets permanently 
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FIG. 4: Left: Dimensionless filament densities along the leading edge, where p c — pj + p c , for example. Right: Total 
dimensionless filament density (p c + pi + p r ) along leading edge. 



capped and no longer extends. Also, no further branching can occur along it. 

(4) Every uncapped filament grows by an additional 100 nm in its initially chosen direction (polymerization). 

(5) Steps (l)-(4) are repeated for each uncapped filament until capped. 

Figure 5 demonstrates output from the simulation for a branching angle of 70 degrees with a = 10°. Note that we 
do not explicitly incorporate a membrane into the simulation and we allow for overlaps of filaments due to the thin, 
third dimension. Also, unless specified otherwise, the time step in the simulation is 0.30 seconds, assuming a constant 
G-actin concentration of 10 /iM, the branch rate is 33.33 s~ 1 iim^ 1 and the capping rate is 0.83 s _1 [Bl l32j. 

Growth: We first investigate the optimal relation between 9 and ip. In keeping with the mean field analysis, we 
compute the average number of uncapped filaments generated each time step, denoted as G, with an upper bound 
of 1000 filaments. Note that here we do not distinguish between the two populations, mother and daughter. If the 
average number of uncapped filaments grows with time, the growth is exponential. Of course, eventually, the system 
reaches a steady state presumably due to a finite amount of Arp2/3 or other mechanisms. To study the approach to 
steady state, one must incorporate recycling of G-actin monomers via depolymerization, severing and debranching as 
well. Such mechanisms have been explored by the autocatalytic model developed by Carlsson [34[ . In the autocatalytic 
model, there is an initial overshoot in the number of filaments that has now been observed experimentally [35j] . 

We present measurements of G, averaged over 4000 samples, as a function of 9. Unless otherwise specified, / = 25 
nanometers. We observe in Figure 6 that the optimal relation, 9* = ±ip/2, holds in the two-dimensional simulations. 
We also see evidence of the subdominant population of filaments for 9 > ip. For ip = 70°, the growth at 9 = 90° 
is more comparable to the growth at 9* than for ip = 80°. Again, perhaps the subdominant pattern for 9 = 90° 




FIG. 5: Discrete simulation output with ip — 70°, 6 = 35°, and a — 10°. The length of the horizontal bar is 1 micron. 
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FIG. 6: Left: Growth rate for different branch angles. Right: Growth rate for for c = 0.67 s 1 (and so more growth than for 
c = 0.83 s- 1 ). 



contributes to increased spreading and/or overlaps between generations following the initial filaments. For a smaller 
capping rate, more growth occurs as evidenced in Figure 6, though the same optimal relation holds, as expected. 
When fluctuations are added to the branch angles, the distinguishing feature of zero growth at 9* = tp remains robust 
as expected. Moreover, G broadens near the maximum. See Figure 7. Broadening was also observed in the mean 
field simulations with noise. 

Overlaps: Urban and collaborators use electron tomography to observe many more overlaps than branches in 
lamellipodia [l3|. Their technique allows one to probe the three-dimensional aspect of the cytoskclcton such that 
filaments that appeared to be branches in two-dimensional electron micrograph images turn out to be overlapping 
filaments. Based on the prevalence of overlaps, Urban and collaborators proposed a new model for the structuring of 
lamellipodia [l3j . Arp2/3 nucleates new filaments near the membrane (just as dimerization does) such that there is 
no pre-existing filament and, hence, no memory of its orientation. In other words, there is no branching. 

However, we would like to point out that the prevalence of overlaps does not rule out a branched model. In fact, 
the existence of overlaps is rather natural in a branched model. If each subsequent generation of branches becomes 
exponentially smaller in length, then there will be no overlaps. This is how one embeds a Bethc lattice-a tree graph — 
in a plane such that there are no overlaps. See Figure 8. If this exponential decrease in length with each generation 
does not occur, then overlaps are expected. While the original electron micrographs indicate that the filament length 
increases further back towards the cell body, we do not expect an exponential increase. Therefore, there will be 
crossings/overlaps between branches. In fact, the overlaps can be reinforced by cross-links thereby increasing the 
temporary rigidity of the network. To test this idea, we measure the number of overlaps and compute the ratio, 
X, of the number of overlaps to the number of branch points for each particular 9. Sec Figure 9. We see that the 
ratio peaks where there is optimal growth. Moreover, as the capping rate decreases, the filaments grow longer also 
allowing for more overlaps. We compare the branched model with a fixed branch angle (plus small fluctuations) to 
a branched model where the branch angle is uniformly random between 1 and 89 degrees. We do this to disrupt the 
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FIG. 7: Growth rate for 70 degree branching angle with noise. 



11 




FIG. 8: Bethe lattice with coordination number three. The circles indicate a repeating pattern. 



inheritance in orientation of the fixed branch angle. So, Arp2/3 merely nucleates a filament off a pre-existing filament 
with no memory of filament orientation. We model the lack of inheritance with the completely random branch angle 
to capture some aspect of the recently proposed unbranched model. 

We observe that the number of overlaps compared to the number of branch points is rather large (exceeding 10) for 
certain initial filament orientations. Indeed, the notion of many overlaps does not rule out the notion of branching. 
A decrease in the capping rate increases Xi as expected, since the filaments are typically longer. Moreover, \ rather 
small for the random branch angle as compared to the fixed branch angle model. The lack of inheritance reduces the 
number of potential overlaps. The reduction in overlaps should also ring true for a completely non- branched Arp2/3 
nucleation model as proposed by Urban and collaborators. In Figure 9, we also plot the overlap for different branch 
angles. Note that for ip = 70°, \ is approximately the same for those filaments originally oriented at = 35° and for 
9 = 90°, the subdominant population of filaments. 

Filament tip spatial distributions: Finally, using the two-dimensional discrete simulation, we compute the 
spatial distribution of filament tips. See Figure 10. From one filament centered at x = y = 0, we observe spreading 
in the x— direction of the dendritic array by several microns. While there is not much difference between the different 
branch angles, for the random branch angle model, the broadening in the x-direction is enhanced. However, that 
broadening is not supported by a large number of overlaps making the network more susceptible to buckling. As for 
the y-direction, the smaller branch angle allows for more forward growth, as expected, however, the overlap ratio is 
also smaller. For the random branch angle, the growth in the y-direction is the largest, but, again, there is not much 
structural support via overlaps. 

Combining the overlap data with the distribution of a;-data we observe that the system is spreading out in the 
x-direction as well as overlapping. The spreading allows for the construction of focal adhesions with which the cells 
temporarily adhere to the surface. The overlaps enhance structural support. Both features are simultaneously possible 
in a branched model via the dendritic nucleation model. In the absence of the branches, the system cross-links, albeit 
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FIG. 9: Left: The ratio of overlaps to branch points as a function of 6 with a — 10° for two fixed branch angle curves. Right: 
The ratio of overlaps to branch points as a function of 9 for a smaller capping rate (more overlaps). Here, a — 0°. 
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FIG. 10: Left: Spatial distribution of filament tips in the ^-direction (horizontal) for ip — 60°, 70°, 80° and the random branch 
model. Right: Spatial distribution of filament tips in the y-direction (vertical) for ip = 60°, 70°, 80° and the random branch 
model. For both plots, a — 10°. 



not as effectively as a branched model, but does not spread out in the x-direction. Moreover, the proliferation of 
branches (G), from a material standpoint, results in the effective strengthening of the material as the meshwork size, 
the distance between overlaps, decreases with an increasing number of branches. If one were to suspend disassembly 
of the network, this gradation can be modelled via a spatially varying elastic modulus. More specifically, if one were 
to model the quasi-two-dimensional lamellipodia as a thin elastic plate with a spatially varying elastic modulus, then 
(1) the buckling instability softens in that the system does not undergo discontinuous mode changes and (2) the 
system becomes more robust against out-of-plane buckling as the elastic modulus increases along the direction of 
axial compression [36j . Therefore, branching accounts for spreading, reinforcements via overlaps, and gradation. Only 
reinforcement is possible in a nonbranching model for a fixed Arp2/3 concentration. 



IV. DISCUSSION 



Invoking a geometric notion for collision-based branching between globular Arp2/3 and linear, actin filaments, we 
constructed a mean field model for the orientation of actin filaments near the leading edge of a crawling cell. To study 
the model, we applied the approach of Maly and Borisy [nil ] , who constructed and studied an initial mean field model 
with a different physical basis than ours. The Maly and Borisy approach [To[ invoked a population biology framework 
with branching corresponding to birth and capping corresponding to death. More specifically, they used the Fisherian 
criterion l23l. [23 for maximal reproduction as an optimization condition on the filament orientation. Similar to Maly 
and Borisy jlOj . we found consistency with previous measurements of the distribution of filament orientation with 
respect to the leading edge. In particular, the two, well-defined peaks in the distribution at 9* = ±35° = ±f/;/2 
coincide with the optimal relation, assuming ip = 70°. The fact that both our kinetic model and Maly/Borisy 
model (lOj obtain the same optimal relation despite the differing kinetic assumptions, even in the presence of noise, 
is interesting and calls for further differentiation. 

Our Arp2/3-actin collision-based model predicts a subdominant population of filaments that may account for 
recent measurements on the distribution of filament orientation, which are in apparent contradiction with the earlier 
measurements (26| . The more recent experiment reported a more broad distribution in filament orientation than 
previously measured. Moreover, the subdominant population of filaments may be invoked to more accurately model 
the filament density along the leading edge. Earlier modelling of the filament density demonstrated that larger filament 
densities in the center required smaller filament densities on the sides of lamellipodia [ill f33| . This requirement is not 
so consistent with observation, however. By extending the earlier filament density model to include the subdominant 
population of filaments, this requirement has been relaxed such that the revised filament density model results are 
more consistent with observations of "excess" filament density at the sides of the leading edge. 

To go beyond mean field and study both the positional and orientational degrees of freedom of the actin network 
in its initial growth phase, we implemented a two-dimensional, kinetic simulation. The mean field optimization 
condition persists in the two-dimensional simulation, at least for small fluctuations in the branch angle. It would 
be interesting to extend our collision-based two-dimensional model to include debranching, dcpolymcrization and 
severing so that we can analyze the approach to steady state and compare our results to the autocatalytic model 
developed by Carlsson J34|, which was recently verified experimentally (3oT |. 
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Very recent observations of lamellipodia in motile cells via electron tomography reported many more overlaps 
between filaments than previously estimated using two-dimensional electron micrograph images [l3j. Urban and 
collaborators fl3| used this observation to dispute the dendritic nucleation mdoel and propose a new model of un- 
branched filament nucleation for lamellipodia construction. However, our measurements of the ratio of overlaps to 
branch points are of order 10 using a branched model. For a branched model to have no overlaps, the filament lengths 
must be exponentially decreasing in length with each generation, i.e. the planar embedding of a Bethe lattice in 
two-dimensions. The rather large ratio of overlaps to branch points actually supports the dendritic nucleation model 
with its inherited branch angle. The inheritance increases the potential for overlaps and the pairing of filaments. In 
fact, the pairing observation was also used by Urban and collaborators [l~3| as a mark against the dendritic nucleation 
model. We must also point out that branching also promotes spreading of the network to more readily assemble focal 
adhesions, and gradation to make the network less susceptible to out-of-plane buckling. 

There exists other evidence for an unbranched model for lamellipodia reconstruction promoting cell motility. Based 
on experimental observation, Brieher and collaborators [37J proposed an initial branching motility phase followed by a 
bundled-actin motility phase, facilitated by facsin or other cross-linking proteins. The notion of a filopodia-dominated 
phase of motility cannot be ruled out and may be one of many phases of cell motility. However, reconstituted 
experiments with Arp2/3-actin-fascin demonstrated that Arp2/3 is excluded from the bundling regions [H|. Recent 
modeling supports this notion [3!| . Unfortunately, Urban and collaborators [l3| were unable to determine the spatial 
location of the Arp2/3 in their experiments. We must also point out that the proposal of unbranched Arp2/3 
nucleation implies that the filament density along the leading edge depends purely on the Arp2/3 and not on the 
pre-existing filament density. If the concentration of Arp2/3 is reasonably uniform along the leading edge, then the 
filament density profile will also be reasonably uniform near the center of the leading edge. Some such profiles were 
observed in "rough" crawling cells [ll| . 

Finally, while we have addressed the optimization between the filament orientation and the branch angle, we have 
not addressed the optimization of the branch angle itself. Why does ip « 70°? From the results of this work, only 
when ip > 60° is the redundant second optimum removed, thus paving the way for a suboptimal orientation whose 
center axis is not 9 = 0°. These off-axis populations allow for further spreading and overlapping. We also observe 
that the ratio of overlaps to branch points is approximately the same for the optimized orientation of 9 = ±35° as 
well as 9 = ±90° so that there is an elastic similarity between the two types of orientations. Other speculations as to 
why tp si 70° may be rooted in structural optimization and the like. 

One of us (JMS) would like to acknowledge the hospitality of the Aspen Center of Physics where part of this work 
was completed and finanical support from DMR-0654373. Both authors would like to acknowledge helpful discussion 
with T. Svitkina. 
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